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Abstract 

An expression is derived for the classical free energy difference between two 
configurations of a system, in terms of an ensemble of finite-time measure- 
ments of the work performed in parametrically switching from one configura- 
tion to the other. Two well-known equilibrium identities emerge as limiting 
cases of this result. 
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Consider a finite classical system in contact with a heat reservoir. A central concept 
in thermodynamics is that of the work performed on such a system, when some external 
parameters of the system are made to change with time. (These parameters may represent, 
for instance, the strength of an external field, or the volume of space within which the system 
is confined, or, more abstractly, some particle-particle interactions which are turned on or off 
during the course of a molecular dynamics simulation.) When the parameters are changed 
infinitely slowly along some path 7 from an initial point A to a final point B in parameter 
space, then the total work W performed on the system is equal to the Helmholtz free energy 
difference AF between the initial and final configurations W = AF = F^ — F^. By 
contrast, when the parameters are switched along 7 at a finite rate, then W will depend on 
the microscopic initial conditions of the system and reservoir, and will on average exceed 
AF: 

W > AF. (1) 

Here and in Eq.|^ below, the overbar denotes an average over an ensemble of measurements 
of W, where each measurement is made after first allowing the system and reservoir to 
equilibrate at temperature T, with the parameters fixed at A. (The path 7 from A to B, 
and the rate at which the parameters are switched along this path, remain unchanged from 
one measurement to the next.) The difference W — AF is just the dissipated work, Wdiss, 
associated with the increase of entropy during an irreversible process. 

Eq.|l] is an inequality. By contrast, the new result derived in this paper is the following 
equality: 

exp -f3W = exp -f3AF, (2a) 

or, equivalently, 

AF = -p-^ In exp -(3W, (2b) 

where f3 = l/ksT. This result, which allows one to extract equilibrium information (the 
free energy difference AF) from the ensemble of non- equilibrium (finite-time) measurements 
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described above, is independent of both the path 7 from A to B, and the rate at which the 
parameters are switched along the path. 

Before proceeding with the proof of Eq.^, we estabhsh notation, and then relate Eq.|^ 
to two well-known equilibrium identities for AF. Since we have fixed our attention on a 
particular path 7 in parameter space, it will be convenient to henceforth view the system 
as parametrized by a single quantity A, which increases from to 1 as we travel from A to 
B along 7. Let z = (q, p) denote a point in the phase space of the system, and let Hx{z) 
denote the Hamiltonian for the system, parametrized by the value of A. Next, let Z\ denote 
the partition function, let (■ ■ ■)x denote a canonical average, and let Fx = —(3~^ In Zx denote 
the free energy, all with respect to the Hamiltonian Hx and the temperature T. We are 
interested in the following scenario, which we will refer to as "the switching process": the 
system evolves, in contact with a heat reservoir, as the value of A is switched from to 
1, over a total switching time tg- Without loss of generality, assume a constant switching 
rate, A = t^^. For a given realization of the switching process, the evolution of the system 
is described by a (stochastic) trajectory z(t), and the work performed on the system is the 
time integral of A dHx/dX along this trajectory: 

Now imagine an ensemble of realizations of the switching process (with 7 and tg fixed), with 
the microscopic initial conditions for the system and reservoir generated from a thermal 
equilibrium ensemble at temperature T. Then W may be computed separately for each 
trajectory z(t) in the ensemble, and the overbars appearing in Eqs.|l] and indicate an 
average over the distribution of values of W thus obtained. 

In the limiting cases of infinitely slow and infinitely fast switching of the external 
parameters, we know explicitly the ensemble distribution of values of W, and thus can 
readily check the validity of our central result. In the slow limit {tg —>■ 00), the system 
is in quasi-static equilibrium with the reservoir throughout the switching process, hence 
W = Jq dX {dHx/dX)x for every trajectory in the ensemble. Eq.|^ then reduces to: 
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In the opposite limit {tg 0), the switching of the Hamiltonian is instantaneous, and so 
the work performed is simply W = Hi — Hq = AH, evaluated at the initial conditions 0. 
Since we have a canonical distribution of initial conditions, Eq.|2B| becomes, in this case: 



AF = -13-^ In (exp -[SAH) . (5) 



These two results, Eqs.H and ^, are well-established identities for the free energy difference 
AF Note that both give AF in terms of equilibrium (canonical) averages. By contrast, 
in the intermediate case of finite t^, our ensemble of trajectories lags behind the equilibrium 
distribution in phase space as Hx changes with time. In this sense, Eq.|^ is an explicitly 
non-equilibrium result. 

To prove our central result, it is instructive to first consider what happens when there is 
no heat reservoir during the switching process. The evolution of the system is then described 
by a deterministic trajectory 7i{t) which evolves under Hx{z), as A changes from to 1 over 
a time tg. Consider an ensemble of such trajectories, defined by a canonical distribution of 
initial conditions (at a temperature T) . This ensemble is described by a phase space density 
/(z,t) which satisfies /(z,0) = Zq^ exp ~i3Hq{z), and which evolves under the Liouville 
equation, df/dt + {f,H\} = 0, with A = A(t) = t/ts. Here, {■,■} denotes the Poisson 
bracket. Since the evolution is deterministic, a particular trajectory in this ensemble is 
uniquely specified by single point: there is exactly one trajectory which passes through a 
given z at time t. This means we can define a "work accumulated" function w{z,t), as 
follows. For the trajectory which passes through the point z at time t, w{z,t) is the work 
performed on that trajectory (the time integral of XdH\/dX) up to time t. Since the total 



work W is just the work accumulated up to time tg (Eq.|D, the ensemble average exp —/3W 
may be expressed as 

exp—/3W = J dz f{z,ts) exp —(3w{z,ts)- (6) 
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Now, the work done on an isolated Hamiltonian system is equal to the change in its energy. 
Thus, w{z,t) = Hx{z) —Ho{zo), where zq = zo(z,t) is the initial condition for the trajectory 
which passes through z at time t, and A = A(t). Furthermore, Liouville's theorem tells 
us that phase space density is conserved along any trajectory, hence /(z,t) = /(zo,0) = 
Zq^ exp —/SHq^zq). Combining these results immediately gives 



Since AF = —{3^^ ln(Zi/Zo), we have established the validity of Eq.^ for the case in which 
the system is isolated during the switching process. 

Now consider the situation in which the system is coupled to a reservoir. We assume that 
the system of interest and the reservoir together constitute a larger, isolated Hamiltonian 
system. Let z' denote a point in the phase space of the reservoir, let 7i(z') be the Hamiltonian 
for the reservoir alone, and let y = (z, z') denote a point in the full phase space of system and 
reservoir. Motion in the full phase space is deterministic, and governed by a Hamiltonian 
= Hxiz) + 7i(z') + hint{z,z'), where the interaction term hint couples the system of 
interest to the reservoir. Let Yx be the partition function for Gx- We explicitly assume the 
reservoir to be large enough, and the interaction energy hint small enough 0, that when A is 
held fixed the system of interest samples its phase space according to the Boltzmann factor 
Q-l3Hx{z)_ Now imagine that, at t = 0, we populate the full phase space with a canonical 
distribution of initial conditions, using the Boltzmann factor e"'^'-^"'-^-'. (This corresponds 
to allowing the coupled system and reservoir to equilibrate at temperature T, before each 
realization of the switching process.) From this ensemble of initial conditions, an ensemble 
of trajectories y(t) evolves deterministically under Gx, as A switches from to 1. Since 
the system of interest and reservoir together constitute an isolated Hamiltonian system, the 
work W performed on the system of interest is equal to the change in the total energy of the 



/(z, t) exp -f]w{z, t) = Zq^ exp -(3Hx{z). 



(7) 



Eq.p then becomes 




(8) 
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system and reservoir: W = Giyy(ts)j — G'o(y(0)j- Therefore, applying the analysis of the 
previous paragraph to the situation considered here, with y, Gx, and Yx replacing z, Hx, 
and Zx, respectively, we get 



The right side of Eq.|^ depends only on the initial and final Hamiltonians Gq and Gi, and on 
the temperature T, which means that the ensemble average exp —f3W is independent of the 
switching time ts (and also of the path from A to i? in parameter space). But we already 
know that exp —j3W = exp —(3AF in the limit ts oo, since W = AF for every member 
of the ensemble, in that limiting case. We therefore conclude that 



Eq.^ which tells us that the ensemble average exp —jSW is independent of both 7 and 
tsi is identically true, given the formulation of the problem. However, in going from Eq.^ 
to Eq.0, we invoke a result from quasi-equilibrium statistical mechanics, which relies on 
the assumption of weak coupling (small hint)- Eq.|, therefore, is valid for sufficiently weak 
coupling between the system of interest and the reservoir. This may be seen more directly 
by writing an explicit expression for the ratio Yi/Yq appearing on the right side of Eq.^: 
only if hint may be neglected does this ratio immediately reduce to Zi/Zq (= exp — /3AF). 

Note that the inequality W > AF (Eq.|l]) follows directly from the equality exp —jSW = 
exp — (Eq.|2a|), by application of the mathematical identity exp x > expx 0. This 
establishes W > AF directly from a microscopic, Hamiltonian basis, rather than by invoking 
the increase of entropy. (In the limit — 0, we have W = {AH)q, and Eq.[^ reduces to the 
Gibbs-Bogoliubov-Feynman bound 0, {AH)q > AF.) 

It is also worthwhile to point out that the right side of Eq.^ may be expanded as a sum 
of cumulants (see Eq.[9] of Ref. [|]): 



exp-(3W = Y,/Yo. 



(9) 



exp -pW = exp -pAF 



(10) 



for all values of ts (and all paths 7). This proves our central result, Eq.^ 



00 




AF = E(-/3) 



(11) 



n=l 
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where uon is the n'th cumulant of the ensemble distribution of values of W . If this distribution 
happens to be Gaussian (as may be expected for sufficiently slow switching), then only the 
first two terms survive, and we have 

AF = TF-/3aV2, (12) 

where cr^ = W^ — W is the ensemble variance of W . The dissipated work Wdiss (= W — ^F) 
is then related to the fiuctuations in W by: Wdiss = /?o"^/2. This is a fiuctuation-dissipation 
relation, and has been obtained within the context of numerical simulations by Hermans . 
(A related result for microcanonical ensembles has been derived by Ott [§].) 

The central result of this paper, Eq.|^, makes a concrete prediction regarding the outcome 
of an ensemble of measurements, which in principle is subject to experimental verification. 
Moreover, this result ought to be valid quite generally, provided the coupling to the reservoir 
is sufficiently weak, and quantal effects may be ignored. In practice, however, the applica- 
bility of Eq.|^ may be severely limited by the following considerations. If the fiuctuations 
in W from one measurement to the next are much larger than ksT (i.e. if a ^ then 
the ensemble average of exp —f3W may be dominated by values of W many standard de- 
viations below W . Since such values of the work represent statistically very rare events, it 



would require an unreasonably large number of measurements of W to determine exp —f3W 
with good accuracy. Therefore, given a specific system of interest, switching path 7, and 
switching time tg, the fiuctuations in the work W must not be much greater than /c^T, if 
we are to have any hope of verifying Eq.^ experimentally. This condition pretty much rules 
out macroscopic systems of interest. In recent years, however, the direct manipulation of 
nanoscale objects — and the measurement of forces thereon — has become feasible. Such 
systems may offer the best chance for experimentally testing the new result of this paper. 

So far, we have implicitly assumed that our system is coupled to a physical heat reservoir. 
It is interesting, however, to discuss this problem within the context of numerical simulations. 
On a computer, a heat reservoir must somehow be "mocked up". One way to accomplish 



this is with a Nose-Hoover (NH) thermostat [|lOl, or some variant thereof. In its simplest 



form, this method replaces the reservoir with a single variable (; motion in the extended 
phase space (z, () is governed by the NH equations 

|g = p/m , p=-V$a-Cp}^ (13) 
C = {K/Ko - l)/r^ (14) 

[We have assumed a kinetic + potential Hamiltonian: H\ = /2m + $A(q)- The index n 
runs over all D degrees of freedom of the system, K = /2m is the total kinetic energy 
of the system, Kq = (3~^D/2 is the thermal average of K, and r is a parameter which acts 
as a relaxation time.] For A fixed, a trajectory z{t) generated by these equations of motion 
samples phase space according to the Boltzmann factor exp —i3Hx{z), provided that the 
evolution is sufficiently chaotic. 

It is interesting to ask, does Eq.^ remain valid if the system evolves under the NH 
equations, rather than under the influence of a physical reservoir? Let us consider an 
ensemble of initial conditions in the extended phase space, described by the density 

/(z, C, 0) = cZ,' exp -/3Qo(z, C), (15) 

where Qx{z,Q = H\{z) + DC/^t'^/213, and c = (Z}r^/27r)^/^ is a normalization factor. (The 
distribution cZ^^ exp —PQ\ is stationary under the NH equations when A is held fixed, and 
may be viewed as the "canonical" distribution in the extended phase space.) Allowing these 
initial conditions to evolve under the NH equations, as A changes from to 1, we obtain an 
ensemble of trajectories described by a time-dependent density /(z, t). As before, the work 
performed on each member of the ensemble is defined to be the time integral of XdHx/dX. 
We now introduce a "work accumulated" function w{z, (, t), analogous to w{z, t) introduced 
earlier. It is straightforward to establish that 

/(z, C, t) = /(zo, Co, 0) exp \d f C{t')dA (16) 

w{z, C, t) = Qxiz, - Qoizo, Co) + P^'D f at')dt', (17) 

Jo 

where (zq, Co) are the initial conditions associated with the trajectory which passes through 
(z,C) at time t, and the integral /q C*^^' is performed along this trajectory. Then, repeating 
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the steps leading to Eq.||, we again get exp —(3W = exp —PAF, where the overbar now 
denotes an average over our ensemble of NH trajectories. Thus, Eq.^ remains valid (given 
the canonical distribution of initial conditions specified by Eq.|l5D when the system is coupled 
to a Nose-Hoover thermostat, as per Eqs.|13| and |14[ This result is identically true: no weak 
coupling assumption is necessary, nor do we need to assume that the evolution is chaotic. 

It may similarly be established that Eq.^ is valid, without additional assumptions, when 
the thermostat is numerically implemented using the Metropolis Monte Carlo algorithm, 
rather than Nose-Hoover dynamics. In that situation, both the system and the Hamiltonian 
evolve by discrete steps, and the work performed is a sum of changes in Hx, evaluated at 
successive locations of the system in phase space ITlI . 



Numerical simulations of this sort are often used to compute free energy differences 
of physical, chemical or biological interest |jl2|. Typically, a number of simulations of slow 
switching from one configuration to another are performed, and the resulting average work is 
used as an upper bound on AF, as per Eq.|^; reversing direction, a lower bound is established 



11|. The central result of the present paper may be useful in this situation: rather than 



taking the straight average of W, one can instead perform the average of exp —f3W, then 
take the logarithm and multiply by — as per Eq.|2b|. In principle this converges to the 
exact value of AF (rather than to an upper or lower bound) as the number of simulations 
tends to infinity. In practice, however, the same note of caution applies here as in the case 
of coupling to a physical heat bath: if the fluctuations in W from one simulation to the 
next are much larger than ksT, then prohibitively many simulations may be necessary to 



determine exp —(3W with the desired accuracy. Thus, we may expect Eq.g to be useful 
in free energy computations, only if a is not much larger than f3~^. Whether or not this 
condition holds for a given system will depend on factors such as the number of degrees of 
freedom, the switching time tg, the switching path 7, and the numerical implementation of 
the heat bath. 

To summarize, the central result of this paper is an equality which gives the free energy 
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difference AF between two configurations A and S of a classical, parameter-dependent 
system, in terms of an ensemble of finite-time measurements of the work performed on the 
system as it is switched from A to B. The derivation of this result relies on the assumption 
of weak coupling between system and reservoir, but otherwise follows directly from the 
properties of Hamilton's equations. Two well-known equilibrium identities for AF, Eqs.^ 
and I, emerge as limiting cases of this more general, non-equilibrium result. Practical 
considerations in all likelihood limit the applicability of Eq.^ to systems of no more than a 
moderate number of degrees of freedom (e.g., nanoscale systems). Finally, the equality may 
be useful when numerical simulations of thermostatted systems are used to compute free 
energy differences. 
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